#!/usr/bin/perl
# 2013/08/12 instead of association frequency, calculate observed / expected
# 2013/05/24 Make graph of distance vs association frequnecy

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;

if(@ARGV != 4 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [map file] -o [output file]\n";
}

my %opt;
getopts("i:o:", \%opt);
my $FILE_map = $opt{i};
my $FILE_out = $opt{o};

my $Resolution = 100;

my %Len_chr = (
	'chr1'	=>	5672121,
	'chr2'	=>	4615470,
	'chr3'	=>	2493767,
);

#---------------------------------------
# collect information
#---------------------------------------
my %observed;
my $total_observed = 0;
my $fh_map = IO::File->new($FILE_map) or die "cannot open $FILE_map: $!";
while($_ = $fh_map->getline()){
	s/\r?\n//;
	my ($id, $chr1, $loc1, $direction1, $mapQ1, $hinID1, $hinLoc1, $chr2, $loc2, $direction2, $mapQ2, $hinID2, $hinLoc2, $type, $num, $locations) = split /\t/;

	if($chr1 ne $chr2){
		next;
	}

	if($type ne 'U:U'){
		next;
	}


	my $distance= abs($loc1 - $loc2);
	my $length = length($distance);

	my $category = int($distance / $Resolution) * $Resolution;

	$observed{$category}++;
	$total_observed++;
}
$fh_map->close();


#---------------------------------------
# possible number calculation
#---------------------------------------
my %possible;
my $total_possible = 0;
for(my $distance = 0; $distance < $Len_chr{'chr1'}; $distance += $Resolution){
	my $n = 0;
	foreach my $chr(qw(chr1 chr2 chr3)){
		if($Len_chr{$chr} - $distance > 0){
			$n += int(($Len_chr{$chr} - $distance) / $Resolution);
		}
	}
	$possible{$distance} = $n;
	$total_possible += $n;
}



#---------------------------------------
# output distance vs freqeuncy
#---------------------------------------
my $fh_out = IO::File->new($FILE_out, 'w') or die "cannot write $FILE_out: $!";
foreach my $distance(sort{$a <=> $b} keys %observed){
	if(exists $possible{$distance}){
		$fh_out->printf("$distance\t$observed{$distance}\t$possible{$distance}\t%.4f\n", $observed{$distance} / $possible{$distance} /($total_observed / $total_possible));
	}else{
		$fh_out->printf("$distance\t$observed{$distance}\t0\t%.4f\n", $observed{$distance} /($total_observed / $total_possible));
	}
}
$fh_out->close();